Data prep

suppressWarnings(library(MendelianRandomization))
suppressWarnings(library(digest))
suppressWarnings(library(tidyr))                         
suppressWarnings(library(MRPRESSO))
suppressWarnings(library(rowr))
suppressWarnings(library(ggplot2))
suppressWarnings(library(dplyr))

setwd('C:/Users/marta/Desktop/IC/Thesis/Genetic Variants')
diabetes_all<-read.csv('C:/Users/marta/Desktop/IC/Thesis/Genetic Variants/Diabetes/Newrf_diabetes_instruments.csv', header = TRUE, stringsAsFactors = FALSE)
#diabetes_all<-diabetes_all[diabetes_all$F.statistic>10,] Use command to remove weak instruments


breast<-readRDS('C:/Users/marta/Desktop/IC/Thesis/Genetic Variants/Breast Cancer/breast.rds')
prostate<-readRDS('C:/Users/marta/Desktop/IC/Thesis/Genetic Variants/Prostate Cancer/prostate.rds')
breast= separate(data=breast, col=phase3_1kg_id, into=c('SNP','rest'), sep= "\\:") #Create SNP column in breast data

#Merge diabetes and breast data
index_b<-match(breast$SNP,diabetes_all$SNP..rs.id.)
diabetes_b<-cbind(diabetes_all[index_b,],breast)

#ALign effect alleles in 2 datasets
data_in_ea = diabetes_b$Effect.Allele
data_in_nea = diabetes_b$Other.Allele

ea_new = diabetes_b$a1
nea_new = diabetes_b$a0


stay=as.numeric(data_in_ea == ea_new)  #stay = 1 if the same
swap=as.numeric(data_in_nea == ea_new) 

#checks
table(stay)
## stay
##  0  1 
## 57 54
table(swap)
## swap
##  0  1 
## 54 57
stay*swap == 0
##   [1] TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE
##  [15] TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE
##  [29] TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE
##  [43] TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE
##  [57] TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE
##  [71] TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE
##  [85] TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE
##  [99] TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE
swap_vec = rep(1, length(swap))
swap_vec = (-2 * swap) + swap_vec

#align the beta vector in the new data
diabetes_b$bcac_onco_icogs_gwas_beta =swap_vec*as.numeric(diabetes_b$bcac_onco_icogs_gwas_beta)
diabetes_b$Effect = swap_vec * as.numeric(diabetes_b$Effect)
####################################################################################

#Merge diabetes and prostate data

index_p<-match(prostate$SNP,diabetes_all$SNP..rs.id.)
diabetes_p<-cbind(diabetes_all[index_p,],prostate)

colnames(diabetes_p)[23]<-'Effect_p'
colnames(diabetes_p)[24]<-'StdErr_p'

#ALign effect alleles in 2 datasets
data_in_ea = diabetes_p$Effect.Allele
data_in_nea = diabetes_p$Other.Allele

ea_new = diabetes_p$Allele1
nea_new = diabetes_p$Allele2


stay=as.numeric(data_in_ea == ea_new)  #stay = 1 if the same
swap=as.numeric(data_in_nea == ea_new) 

#checks
table(stay)
## stay
##  0 
## 88
table(swap)
## swap
##  0 
## 88
stay*swap == 0
##  [1] TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE
## [15] TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE
## [29] TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE
## [43] TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE
## [57] TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE
## [71] TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE
## [85] TRUE TRUE TRUE TRUE
swap_vec = rep(1, length(swap))
swap_vec = (-2 * swap) + swap_vec

#align the beta vector in the new data
diabetes_p$Effect_p =swap_vec*as.numeric(diabetes_p$Effect_p)
diabetes_p$Effect = swap_vec * as.numeric(diabetes_p$Effect)
# remove NAs, we believe that they are missing at random
diabetes_b<-na.omit(diabetes_b)
diabetes_p<-na.omit(diabetes_p)

#Get rid of duplicated columns
diabetes_b <- diabetes_b[, !duplicated(colnames(diabetes_b))]
diabetes_p <- diabetes_p[, !duplicated(colnames(diabetes_p))]
#Load data of rs ids of grouped SNPs from PhenoScanner query

#Load data of rs ids of grouped SNPs of 3 molecular pathways
Lipids1<-read.csv('C:/Users/marta/Desktop/IC/Thesis/Lipids1.csv', stringsAsFactors = FALSE, header=FALSE)
Insulin1<-read.csv('C:/Users/marta/Desktop/IC/Thesis/Insulin1.csv', stringsAsFactors = FALSE, header=FALSE)
Glucose1<-read.csv('C:/Users/marta/Desktop/IC/Thesis/Glucose1.csv', stringsAsFactors = FALSE, header=FALSE)

MR breast cancer

MRInputObject_breast <- mr_input(bx = diabetes_b$Effect,
                          bxse = diabetes_b$StdErr,
                          by = diabetes_b$bcac_onco_icogs_gwas_beta,
                          byse = diabetes_b$bcac_onco_icogs_gwas_se, outcome = 'breast cancer', exposure='diabetes',snps=diabetes_b$SNP)

MRAllObject_all_breast <- mr_allmethods(MRInputObject_breast, method = "all")

Breast cancer MR effect estimates:

MRAllObject_all_breast
##                     Method Estimate Std Error 95% CI        P-value
##              Simple median   -0.004     0.019  -0.041 0.033   0.835
##            Weighted median    0.011     0.018  -0.024 0.045   0.539
##  Penalized weighted median   -0.004     0.018  -0.039 0.031   0.825
##                                                                    
##                        IVW    0.000     0.017  -0.033 0.033   0.998
##              Penalized IVW   -0.007     0.013  -0.033 0.019   0.581
##                 Robust IVW    0.005     0.023  -0.040 0.050   0.819
##       Penalized robust IVW   -0.008     0.016  -0.040 0.024   0.612
##                                                                    
##                   MR-Egger    0.029     0.035  -0.040 0.098   0.407
##                (intercept)   -0.003     0.003  -0.009 0.003   0.344
##         Penalized MR-Egger    0.003     0.031  -0.058 0.064   0.926
##                (intercept)   -0.001     0.002  -0.006 0.004   0.743
##            Robust MR-Egger    0.046     0.066  -0.083 0.176   0.484
##                (intercept)   -0.004     0.005  -0.013 0.006   0.447
##  Penalized robust MR-Egger   -0.005     0.041  -0.085 0.076   0.908
##                (intercept)    0.000     0.003  -0.006 0.006   0.939
mr_plot(MRInputObject_breast,
        error = TRUE, orientate = FALSE, line = "ivw")
mr_plot(MRAllObject_all_breast)

MR-PRESSO

#MR-PRESSO breast

input = data.frame(betaY=diabetes_b$bcac_onco_icogs_gwas_beta, seY=diabetes_b$bcac_onco_icogs_gwas_se,beta_diabetes=diabetes_b$Effect, se_diabetes=diabetes_b$StdErr,row.names= NULL)

mrpresso_out = mr_presso(BetaOutcome = "betaY", BetaExposure = c("beta_diabetes"), SdOutcome = "seY", SdExposure="se_diabetes", data=input, OUTLIERtest = TRUE, DISTORTIONtest = TRUE, NbDistribution = 1000)
## Warning in mr_presso(BetaOutcome = "betaY", BetaExposure =
## c("beta_diabetes"), : Outlier test unstable. The significance threshold
## of 0.05 for the outlier test is not achievable with only 1000 to
## compute the null distribution. The current precision is <0.111. Increase
## NbDistribution.
mrpresso_out$`MR-PRESSO results`$`Global Test`$`RSSobs`
## [1] 459.3003
mrpresso_out$`MR-PRESSO results`$`Global Test`$Pvalue
## [1] "<0.001"
mrpresso_out$`Main MR results`
##        Exposure       MR Analysis Causal Estimate         Sd       T-stat
## 1 beta_diabetes               Raw    4.022679e-05 0.01704120  0.002360562
## 2 beta_diabetes Outlier-corrected   -7.200820e-03 0.01400442 -0.514181786
##     P-value
## 1 0.9981208
## 2 0.6082263

MR-PRESSO Distortion test

mrpresso_out$`MR-PRESSO results`$`Distortion Test`$`Distortion Coefficient`
## beta_diabetes 
##      100.5586
mrpresso_out$`MR-PRESSO results`$`Distortion Test`$Pvalue
## [1] 0.228

MR-PRESSO Outliers detected

out_indices = mrpresso_out$`MR-PRESSO results`$`Distortion Test`$`Outliers Indices`
outliers = diabetes_b$SNP[out_indices]

# creating a scatterplot
ggplot(data = diabetes_b, aes(x = diabetes_b$Effect, y = diabetes_b$bcac_onco_icogs_gwas_beta, color =ifelse(diabetes_b$SNP%in%outliers,yes='red',no='blue'))) + 
  geom_point(size=1.5) +
  labs(title = "", x = "Genetic associations with T2D", y = "Genetic associations with breast cancer", color = "Genetic variants\n") +
  scale_color_manual(labels = c("Rest of SNPs", "outliers"), values = c("blue", "red")) +
  theme_bw() +
  theme(axis.text.x = element_text(size = 10), axis.title.x = element_text(size = 12),
        axis.text.y = element_text(size = 10), axis.title.y = element_text(size = 12),
        plot.title = element_text(size = 20, face = "bold", color = "darkgreen"))

Stratified MR analyses

Insulin

Scatter plot of grouped SNPs

# creating a scatterplot
ggplot(data = diabetes_b, aes(x = diabetes_b$Effect, y = diabetes_b$bcac_onco_icogs_gwas_beta, color =ifelse(diabetes_b$SNP%in%Insulin1$V1,yes='red',no='blue'))) + 
  geom_point(size=2) +
  labs(title = "MR analysis\n", x = "Genetic associations with T2D", y = "Genetic associations with breast cancer", color = "Genetic variants\n") +
  scale_color_manual(labels = c("Rest of SNPs", "Insulin SNPs"), values = c("blue", "red")) +
  theme_bw() +
  theme(axis.text.x = element_text(size = 14), axis.title.x = element_text(size = 16),
        axis.text.y = element_text(size = 14), axis.title.y = element_text(size = 16),
        plot.title = element_text(size = 20, face = "bold", color = "darkgreen"))

index<-match(diabetes_b$SNP,Insulin1$V1)
data_b_insulin<-na.omit(diabetes_b[index,])
MRInputObject_breast_insulin <- mr_input(bx = data_b_insulin$Effect,
                          bxse = data_b_insulin$StdErr,
                          by = data_b_insulin$bcac_onco_icogs_gwas_beta,
                          byse = data_b_insulin$bcac_onco_icogs_gwas_se, outcome = 'breast cancer', exposure='diabetes',snps=data_b_insulin$SNP)

MRAllObject_all_breast_insulin <- mr_allmethods(MRInputObject_breast_insulin, method = "all")
MRAllObject_all_breast_insulin
##                     Method Estimate Std Error 95% CI         P-value
##              Simple median   -0.073     0.121  -0.311  0.164   0.544
##            Weighted median   -0.053     0.080  -0.209  0.103   0.507
##  Penalized weighted median   -0.154     0.075  -0.301 -0.007   0.040
##                                                                     
##                        IVW   -0.010     0.085  -0.178  0.157   0.903
##              Penalized IVW   -0.108     0.074  -0.254  0.038   0.148
##                 Robust IVW   -0.017     0.138  -0.288  0.253   0.901
##       Penalized robust IVW   -0.123     0.082  -0.284  0.039   0.136
##                                                                     
##                   MR-Egger    0.012     0.141  -0.265  0.289   0.932
##                (intercept)   -0.002     0.009  -0.020  0.016   0.835
##         Penalized MR-Egger   -0.037     0.134  -0.298  0.225   0.784
##                (intercept)    0.000     0.008  -0.016  0.016   0.993
##            Robust MR-Egger    0.008     0.151  -0.287  0.303   0.959
##                (intercept)   -0.002     0.007  -0.015  0.011   0.782
##  Penalized robust MR-Egger   -0.039     0.121  -0.276  0.197   0.745
##                (intercept)    0.000     0.006  -0.012  0.012   0.982
mr_plot(MRInputObject_breast_insulin,
        error = TRUE, orientate = FALSE, line = "ivw")
mr_plot(MRAllObject_all_breast_insulin)

Glucose

Scatter plot of grouped SNPs

# creating a scatterplot
ggplot(data = diabetes_b, aes(x = diabetes_b$Effect, y = diabetes_b$bcac_onco_icogs_gwas_beta, color =ifelse(diabetes_b$SNP%in%Glucose1$V1,yes='red',no='blue'))) + 
  geom_point(size=2) +
  labs(title = "MR analysis\n", x = "Genetic associations with T2D", y = "Genetic associations with breast cancer", color = "Genetic variants\n") +
  scale_color_manual(labels = c("Rest of SNPs", "Glucose SNPs"), values = c("blue", "red")) +
  theme_bw() +
  theme(axis.text.x = element_text(size = 14), axis.title.x = element_text(size = 16),
        axis.text.y = element_text(size = 14), axis.title.y = element_text(size = 16),
        plot.title = element_text(size = 20, face = "bold", color = "darkgreen"))

index<-match(diabetes_b$SNP,Glucose1$V1)
data_b_glucose<-na.omit(diabetes_b[index,])
MRInputObject_breast_glucose <- mr_input(bx = data_b_glucose$Effect,
                          bxse = data_b_glucose$StdErr,
                          by = data_b_glucose$bcac_onco_icogs_gwas_beta,
                          byse = data_b_glucose$bcac_onco_icogs_gwas_se, outcome = 'breast cancer', exposure='diabetes',snps=data_b_glucose$SNP)

MRAllObject_all_breast_glucose <- mr_allmethods(MRInputObject_breast_glucose, method = "all")
MRAllObject_all_breast_glucose
##                     Method Estimate Std Error 95% CI        P-value
##              Simple median   -0.005     0.118  -0.237 0.227   0.965
##            Weighted median    0.103     0.068  -0.031 0.237   0.131
##  Penalized weighted median    0.209     0.068   0.075 0.343   0.002
##                                                                    
##                        IVW    0.041     0.090  -0.135 0.217   0.650
##              Penalized IVW    0.009     0.064  -0.117 0.135   0.883
##                 Robust IVW    0.044     0.101  -0.155 0.242   0.665
##       Penalized robust IVW    0.014     0.103  -0.189 0.216   0.895
##                                                                    
##                   MR-Egger    0.094     0.159  -0.217 0.405   0.552
##                (intercept)   -0.004     0.010  -0.023 0.015   0.676
##         Penalized MR-Egger    0.073     0.100  -0.124 0.269   0.469
##                (intercept)    0.000     0.007  -0.013 0.013   0.985
##            Robust MR-Egger    0.095     0.142  -0.184 0.373   0.506
##                (intercept)   -0.004     0.009  -0.021 0.013   0.651
##  Penalized robust MR-Egger    0.094     0.178  -0.255 0.443   0.597
##                (intercept)   -0.001     0.008  -0.017 0.015   0.916
mr_plot(MRInputObject_breast_glucose,
        error = TRUE, orientate = FALSE, line = "ivw")
mr_plot(MRAllObject_all_breast_glucose)

Lipids

index<-match(diabetes_b$SNP,Lipids1$V1)
data_b_lipid<-na.omit(diabetes_b[index,])

Scatter plot of grouped SNPs

# creating a scatterplot
ggplot(data = diabetes_b, aes(x = diabetes_b$Effect, y = diabetes_b$bcac_onco_icogs_gwas_beta, color =ifelse(diabetes_b$SNP%in%Lipids1$V1,yes='red',no='blue'))) + 
  geom_point(size=2) +
  labs(title = "MR analysis\n", x = "Genetic associations with T2D", y = "Genetic associations with breast cancer", color = "Genetic variants\n") +
  scale_color_manual(labels = c("Rest of SNPs", "Lipid SNPs"), values = c("blue", "red")) +
  theme_bw() +
  theme(axis.text.x = element_text(size = 14), axis.title.x = element_text(size = 16),
        axis.text.y = element_text(size = 14), axis.title.y = element_text(size = 16),
        plot.title = element_text(size = 20, face = "bold", color = "darkgreen"))

MRInputObject_breast_lipid <- mr_input(bx = data_b_lipid$Effect,
                          bxse = data_b_lipid$StdErr,
                          by = data_b_lipid$bcac_onco_icogs_gwas_beta,
                          byse = data_b_lipid$bcac_onco_icogs_gwas_se, outcome = 'breast cancer', exposure='diabetes',snps=data_b_lipid$SNP)

MRAllObject_all_breast_lipid <- mr_allmethods(MRInputObject_breast_lipid, method = "all")
MRAllObject_all_breast_lipid
##                     Method Estimate Std Error 95% CI        P-value
##              Simple median    0.105     0.089  -0.070 0.280   0.239
##            Weighted median    0.216     0.079   0.062 0.370   0.006
##  Penalized weighted median    0.226     0.082   0.065 0.388   0.006
##                                                                    
##                        IVW    0.111     0.102  -0.089 0.311   0.276
##              Penalized IVW    0.162     0.075   0.016 0.309   0.030
##                 Robust IVW    0.118     0.125  -0.128 0.363   0.348
##       Penalized robust IVW    0.166     0.056   0.056 0.276   0.003
##                                                                    
##                   MR-Egger   -0.003     0.230  -0.453 0.447   0.989
##                (intercept)    0.010     0.018  -0.026 0.047   0.571
##         Penalized MR-Egger    0.093     0.222  -0.342 0.528   0.675
##                (intercept)    0.006     0.017  -0.027 0.039   0.719
##            Robust MR-Egger    0.000     0.243  -0.477 0.478   0.999
##                (intercept)    0.010     0.020  -0.028 0.049   0.597
##  Penalized robust MR-Egger    0.096     0.207  -0.309 0.502   0.641
##                (intercept)    0.006     0.019  -0.030 0.042   0.745
mr_plot(MRInputObject_breast_lipid,
        error = TRUE, orientate = FALSE, line = "ivw")
mr_plot(MRAllObject_all_breast_lipid)

No BMI

BMI=read.csv('C:/Users/marta/Desktop/IC/Thesis/Genetic variants/BMI.csv', stringsAsFactors = FALSE, header=FALSE)
index<-match(diabetes_b$SNP,BMI$V1)
data_b_bmi<-cbind(diabetes_b,BMI[index,])
data_b_bmi$V2[is.na(data_b_bmi$V2)] <- 'NO BMI'
data_nobmi<-data_b_bmi[data_b_bmi$V2=='NO BMI',]

MRInputObject_breast_nobmi <- mr_input(bx = data_nobmi$Effect,
                          bxse = data_nobmi$StdErr,
                          by = data_nobmi$bcac_onco_icogs_gwas_beta,
                          byse = data_nobmi$bcac_onco_icogs_gwas_se, outcome = 'breast cancer', exposure='diabetes',snps=data_nobmi$SNP)

MRAllObject_all_breast_nobmi <- mr_allmethods(MRInputObject_breast_nobmi, method = "all")
MRAllObject_all_breast_nobmi
##                     Method Estimate Std Error 95% CI        P-value
##              Simple median    0.036     0.030  -0.022 0.095   0.221
##            Weighted median   -0.005     0.027  -0.059 0.049   0.856
##  Penalized weighted median   -0.004     0.028  -0.058 0.050   0.871
##                                                                    
##                        IVW   -0.012     0.022  -0.056 0.031   0.583
##              Penalized IVW    0.002     0.021  -0.038 0.043   0.906
##                 Robust IVW   -0.010     0.022  -0.052 0.033   0.649
##       Penalized robust IVW    0.003     0.020  -0.036 0.041   0.898
##                                                                    
##                   MR-Egger   -0.040     0.048  -0.134 0.055   0.408
##                (intercept)    0.002     0.003  -0.004 0.009   0.517
##         Penalized MR-Egger   -0.043     0.046  -0.133 0.048   0.358
##                (intercept)    0.003     0.003  -0.004 0.009   0.374
##            Robust MR-Egger   -0.047     0.041  -0.127 0.034   0.254
##                (intercept)    0.003     0.004  -0.004 0.011   0.400
##  Penalized robust MR-Egger   -0.049     0.039  -0.126 0.029   0.217
##                (intercept)    0.004     0.003  -0.003 0.010   0.296
mr_plot(MRInputObject_breast_nobmi,
        error = TRUE, orientate = FALSE, line = "ivw")
mr_plot(MRAllObject_all_breast_nobmi)

Prostate cancer MR effect estimates:

MRInputObject_prostate <- mr_input(bx = diabetes_p$Effect,
                          bxse = diabetes_p$StdErr,
                          by = diabetes_p$Effect_p,
                          byse = diabetes_p$StdErr_p, outcome = 'prostate cancer', exposure='diabetes',snps=diabetes_p$SNP)

MRAllObject_all_prostate <- mr_allmethods(MRInputObject_prostate, method = "all")
MRAllObject_all_prostate 
##                     Method Estimate Std Error 95% CI        P-value
##              Simple median    0.011     0.024  -0.035 0.058   0.632
##            Weighted median    0.001     0.022  -0.042 0.045   0.947
##  Penalized weighted median   -0.005     0.022  -0.049 0.039   0.819
##                                                                    
##                        IVW   -0.041     0.042  -0.124 0.042   0.333
##              Penalized IVW   -0.010     0.016  -0.042 0.022   0.543
##                 Robust IVW   -0.007     0.018  -0.042 0.028   0.690
##       Penalized robust IVW   -0.012     0.016  -0.043 0.020   0.469
##                                                                    
##                   MR-Egger   -0.113     0.100  -0.309 0.084   0.260
##                (intercept)    0.006     0.007  -0.009 0.021   0.428
##         Penalized MR-Egger   -0.048     0.038  -0.122 0.027   0.213
##                (intercept)    0.003     0.003  -0.003 0.009   0.297
##            Robust MR-Egger   -0.020     0.039  -0.097 0.056   0.603
##                (intercept)    0.001     0.003  -0.004 0.006   0.690
##  Penalized robust MR-Egger   -0.045     0.032  -0.108 0.017   0.155
##                (intercept)    0.003     0.002  -0.002 0.007   0.274
mr_plot(MRInputObject_prostate,
        error = TRUE, orientate = FALSE, line = "ivw")
mr_plot(MRAllObject_all_prostate)

MR-PRESSO prostate

#MR-PRESSO prostate

input = data.frame(betaY=diabetes_p$Effect_p,seY=diabetes_p$StdErr_p,beta_diabetes=diabetes_p$Effect,se_diabetes=diabetes_p$StdErr, row.names= NULL)

mrpresso_out = mr_presso(BetaOutcome = "betaY", BetaExposure = "beta_diabetes", SdOutcome = "seY", SdExposure="se_diabetes", data=input, OUTLIERtest = TRUE, DISTORTIONtest = TRUE, NbDistribution = 1000)
## Warning in mr_presso(BetaOutcome = "betaY", BetaExposure =
## "beta_diabetes", : Outlier test unstable. The significance threshold of
## 0.05 for the outlier test is not achievable with only 1000 to compute
## the null distribution. The current precision is <0.088. Increase
## NbDistribution.
mrpresso_out$`Main MR results`
##        Exposure       MR Analysis Causal Estimate         Sd     T-stat
## 1 beta_diabetes               Raw     -0.04085910 0.04223478 -0.9674276
## 2 beta_diabetes Outlier-corrected     -0.01429709 0.01722820 -0.8298654
##     P-value
## 1 0.3360120
## 2 0.4089954

MR-PRESSO Distortion test

mrpresso_out$`MR-PRESSO results`$`Distortion Test`$`Distortion Coefficient`
## beta_diabetes 
##     -185.7861
mrpresso_out$`MR-PRESSO results`$`Distortion Test`$Pvalue
## [1] 0.044

MR-PRESSO Outliers detected

out_indices = mrpresso_out$`MR-PRESSO results`$`Distortion Test`$`Outliers Indices`
outliers=diabetes_p$SNP[out_indices]


ggplot(data = diabetes_p, aes(x = diabetes_p$Effect, y = diabetes_p$Effect_p, color =ifelse(diabetes_p$SNP%in%outliers,yes='red',no='blue'))) + 
  geom_point(size=1.5) +
  labs(title='',x = "Genetic association with T2D", y = "Genetic association with prostate cancer", color = "Genetic variants\n") +
  scale_color_manual(labels = c("Rest of SNPs", "outliers"), values = c("blue", "red")) +
  theme_bw() +
  theme(axis.text.x = element_text(size = 10), axis.title.x = element_text(size = 12),
        axis.text.y = element_text(size = 10), axis.title.y = element_text(size = 12),
        plot.title = element_text(size = 20, face = "bold", color = "darkgreen"))

MR stratified analysis

Insulin

index<-match(diabetes_p$SNP,Insulin1$V1)
data_p_insulin<-na.omit(diabetes_p[index,])
ggplot(data = diabetes_p, aes(x = diabetes_p$Effect, y = diabetes_p$Effect_p, color =ifelse(diabetes_p$SNP%in%Insulin1$V1,yes='red',no='blue'))) + 
  geom_point(size=2) +
  labs(title = "MR analysis\n", x = "T2D", y = "Prostate cancer", color = "Genetic variants\n") +
  scale_color_manual(labels = c("Rest of SNPs", "Insulin SNPs"), values = c("blue", "red")) +
  theme_bw() +
  theme(axis.text.x = element_text(size = 14), axis.title.x = element_text(size = 16),
        axis.text.y = element_text(size = 14), axis.title.y = element_text(size = 16),
        plot.title = element_text(size = 20, face = "bold", color = "darkgreen"))

MRInputObject_prostate_insulin <- mr_input(bx = data_p_insulin$Effect,
                          bxse = data_p_insulin$StdErr,
                          by = data_p_insulin$Effect_p,
                          byse = data_p_insulin$StdErr_p, outcome = 'prostate cancer', exposure='T2D',snps=data_p_insulin$SNP)

MRAllObject_all_prostate <- mr_allmethods(MRInputObject_prostate_insulin, method = "all")
MRAllObject_all_prostate 
##                     Method Estimate Std Error 95% CI        P-value
##              Simple median    0.060     0.114  -0.164 0.285   0.598
##            Weighted median    0.061     0.095  -0.125 0.247   0.520
##  Penalized weighted median    0.085     0.097  -0.105 0.274   0.383
##                                                                    
##                        IVW    0.008     0.087  -0.163 0.180   0.925
##              Penalized IVW    0.062     0.082  -0.098 0.222   0.447
##                 Robust IVW    0.023     0.124  -0.220 0.267   0.850
##       Penalized robust IVW    0.067     0.080  -0.090 0.223   0.403
##                                                                    
##                   MR-Egger   -0.001     0.178  -0.349 0.347   0.995
##                (intercept)    0.001     0.011  -0.021 0.023   0.950
##         Penalized MR-Egger   -0.001     0.178  -0.349 0.347   0.995
##                (intercept)    0.001     0.011  -0.021 0.023   0.950
##            Robust MR-Egger    0.044     0.331  -0.604 0.693   0.894
##                (intercept)   -0.001     0.012  -0.025 0.023   0.932
##  Penalized robust MR-Egger    0.044     0.331  -0.604 0.693   0.894
##                (intercept)   -0.001     0.012  -0.025 0.023   0.932
mr_plot(MRInputObject_prostate_insulin,
        error = TRUE, orientate = FALSE, line = "ivw")
mr_plot(MRAllObject_all_prostate)

Glucose

index<-match(diabetes_p$SNP,Glucose1$V1)
data_p_glucose<-na.omit(diabetes_p[index,])
ggplot(data = diabetes_p, aes(x = diabetes_p$Effect, y = diabetes_p$Effect_p, color =ifelse(diabetes_p$SNP%in%Glucose1$V1,yes='red',no='blue'))) + 
  geom_point(size=2) +
  labs(title = "MR analysis\n", x = "T2D", y = "Prostate cancer", color = "Genetic variants\n") +
  scale_color_manual(labels = c("Rest of SNPs", "Glucose SNPs"), values = c("blue", "red")) +
  theme_bw() +
  theme(axis.text.x = element_text(size = 14), axis.title.x = element_text(size = 16),
        axis.text.y = element_text(size = 14), axis.title.y = element_text(size = 16),
        plot.title = element_text(size = 20, face = "bold", color = "darkgreen"))

MRInputObject_prostate_glucose <- mr_input(bx = data_p_glucose$Effect,
                          bxse = data_p_glucose$StdErr,
                          by = data_p_glucose$Effect_p,
                          byse = data_p_glucose$StdErr_p, outcome = 'prostate cancer', exposure='diabetes',snps=data_p_glucose$SNP)

MRAllObject_all_prostate <- mr_allmethods(MRInputObject_prostate_glucose, method = "all")
MRAllObject_all_prostate 
##                     Method Estimate Std Error 95% CI        P-value
##              Simple median    0.163     0.110  -0.053 0.379   0.139
##            Weighted median    0.094     0.088  -0.078 0.266   0.285
##  Penalized weighted median    0.098     0.094  -0.086 0.281   0.297
##                                                                    
##                        IVW    0.108     0.086  -0.061 0.276   0.210
##              Penalized IVW    0.216     0.060   0.098 0.334   0.000
##                 Robust IVW    0.137     0.146  -0.150 0.424   0.350
##       Penalized robust IVW    0.211     0.092   0.032 0.391   0.021
##                                                                    
##                   MR-Egger    0.134     0.166  -0.191 0.459   0.418
##                (intercept)   -0.002     0.012  -0.025 0.021   0.849
##         Penalized MR-Egger    0.335     0.109   0.121 0.549   0.002
##                (intercept)   -0.010     0.007  -0.023 0.004   0.164
##            Robust MR-Egger    0.385     0.329  -0.260 1.030   0.242
##                (intercept)   -0.013     0.024  -0.060 0.035   0.599
##  Penalized robust MR-Egger    0.353     0.299  -0.234 0.940   0.238
##                (intercept)   -0.011     0.021  -0.051 0.029   0.587
mr_plot(MRInputObject_prostate_glucose,
        error = TRUE, orientate = FALSE, line = "ivw")
mr_plot(MRAllObject_all_prostate)

Lipids

index<-match(diabetes_p$SNP,Lipids1$V1)
data_p_lipid<-na.omit(diabetes_p[index,])
ggplot(data = diabetes_p, aes(x = diabetes_p$Effect, y = diabetes_p$Effect_p, color =ifelse(diabetes_p$SNP%in%Lipids1$V1,yes='red',no='blue'))) + 
  geom_point(size=2) +
  labs(title = "MR analysis\n", x = "T2D", y = "Prostate cancer", color = "Genetic variants\n") +
  scale_color_manual(labels = c("Rest of SNPs", "Lipid SNPs"), values = c("blue", "red")) +
  theme_bw() +
  theme(axis.text.x = element_text(size = 14), axis.title.x = element_text(size = 16),
        axis.text.y = element_text(size = 14), axis.title.y = element_text(size = 16),
        plot.title = element_text(size = 20, face = "bold", color = "darkgreen"))

MRInputObject_prostate_lipid <- mr_input(bx = data_p_lipid$Effect,
                          bxse = data_p_lipid$StdErr,
                          by = data_p_lipid$Effect_p,
                          byse = data_p_lipid$StdErr_p, outcome = 'prostate cancer', exposure='diabetes',snps=data_p_lipid$SNP)

MRAllObject_all_prostate <- mr_allmethods(MRInputObject_prostate_lipid, method = "all")
MRAllObject_all_prostate 
##                     Method Estimate Std Error 95% CI         P-value
##              Simple median    0.090     0.100  -0.105  0.286   0.366
##            Weighted median    0.089     0.089  -0.085  0.263   0.317
##  Penalized weighted median    0.089     0.095  -0.098  0.276   0.351
##                                                                     
##                        IVW    0.105     0.100  -0.090  0.301   0.291
##              Penalized IVW    0.156     0.079   0.001  0.310   0.048
##                 Robust IVW    0.116     0.133  -0.145  0.377   0.384
##       Penalized robust IVW    0.155     0.068   0.021  0.289   0.023
##                                                                     
##                   MR-Egger    0.204     0.248  -0.282  0.689   0.410
##                (intercept)   -0.009     0.019  -0.046  0.029   0.659
##         Penalized MR-Egger    0.542     0.143   0.262  0.823   0.000
##                (intercept)   -0.026     0.010  -0.046 -0.006   0.011
##            Robust MR-Egger    0.571     0.037   0.498  0.645   0.000
##                (intercept)   -0.025     0.002  -0.029 -0.020   0.000
##  Penalized robust MR-Egger    0.571     0.037   0.498  0.645   0.000
##                (intercept)   -0.025     0.002  -0.029 -0.020   0.000
mr_plot(MRInputObject_prostate_lipid,
        error = TRUE, orientate = FALSE, line = "ivw")
mr_plot(MRAllObject_all_prostate)